Velocity Distribution and the Effect of Wall Roughness in Granular 

Poiseuille Flow 



K. C. Vijayakumar and Meheboob Alam* 
Engineering Mechanics Unit, Jawaharlal Nehru Center for 
Advanced Scientific Research, Jakkur PO, Bangalore 560064, India 
(Dated: February 4, 2008) 

From event-driven simulations of a gravity-driven channel flow of inelastic hard-disks, 
we show that the velocity distribution function remains close to a Gaussian for a wide range 
densities (even when the Knudsen number is of order one) if the walls are smooth and 
the particle collisions are nearly elastic. For dense flows, a transition from a Gaussian to 
a power-law distribution for the high velocity tails occurs with increasing dissipation in 
the center of the channel, irrespective of wall-roughness. For a rough wall, the near-wall 
distribution functions are distinctly different from those in the bulk even in the quasielastic 
limit. 

PACS numbers: 45.70.-n, 47.57.Gc 

Granular materials, a collection of macroscopic particles, are important in many chemical and 
pharmaceutical industries as well as in geophysical contexts (avalanche, sand dune, etc.). In the 
rapid flow regime [1], the theory for flowing granular materials is largely based on the dense gas 
kinetic theory that incorporates the inelastic nature of particle collisions. At the heart of such gas- 
or liquid-state continuum theories lies the concept of 'coarse-graining' over distribution functions 
while making a transition from the particle-level properties to the macro- scale fields. Unlike the 
molecular fluid for which the Maxwell-Boltzmann (Gaussian/Maxwellian) distribution plays the 
role of the 'equilibrium' distribution function, however, the granular fluid does not possess any 
'equilibrium' state yj,|2|] due to the microscopic dissipation of particle collisions. However, there 
are 'non-equilibrium' (driven) steady-states for various canonical granular flow configurations 
for which the Gaussian distribution is the leading-order velocity distribution B2D in appropriate 
limits. A systematic study of distribution functions is, therefore, of interest from the viewpoint of 
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developing constitutive models for granular flows as well as to pinpoint the range of validity of 
any adopted theory. Another important issue that needs attention is the derivation of continuum 
boundary conditions for granular flows yfl where it is generally assumed that the distribution 
function in the near wall-region is the same as that in the bulk which is unlikely to hold as we shall 
show here. 

In 'driven' granular flows, the deviation of velocity distribution from a Gaussian has been stud- 
ied through theory [Q, simulation [fj] and experiment [Q, si Is]. In Ref. [5] it has been 
shown that the velocity fluctuations in a vibrated bed of particles follow a Gaussian distribution 
in the solid phase and a power-law distribution (with an exponent —3) in the fluidized phase. 
For the plane shear flow [2], the velocity distribution function is well fitted by an exponent of 
a second-order polynomial in the norm of the fluctuating velocities with angle-dependent coeffi- 
cients. Theoretical works Q4£| for a randomly heated granular gas, based on the Boltzmann-Enskog 
equation, have predicted velocity distribution functions of the form P(v) ~ exp (—717°), with the 
exponent of high- velocity tails being a = 3/2 which also depends on the level of inelastic dissi- 
pation. Experiments [7] for a granular gas confined between two vertical plates and driven into a 
steady state via vertical vibrations have shown that a ~ 1.55 ± 0.1, for a wide range of frequency 
and amplitude of vibrations. Some recent experiments [8], however, showed that the high- velocity 
tails cannot be described by a 'single' universal exponent. 

In this paper we consider the 'granular' Poiseuille flow which is the gravity-driven flow of gran- 
ular materials through a two-dimensional channel Jlp| . focusing on the 'rapid' flow regime 1 1]. 
The simulated system is a channel of length L along the periodic x-direction and bounded by two 
plane solid walls, parallel to the x-direction, with a separation of width W (along the y-direction). 
The granular material, consisting of N identical rigid and smooth disks of equal mass and diameter 
d, is driven by gravity along the x-direction. The interactions that are allowed are instantaneous 
'dissipative' collisions between pairs of particles and between a particle and the walls, via an 
event-driven algorithm flljj]. The dissipative nature of particle collisions is characterized by the 
coefficient of normal restitution, e n , which is the ratio between the pre- and post-collisional rela- 
tive velocities of the colliding particles. There is no relative tangential velocity since the particles 
are assumed to be smooth. The solid walls are modeled as frictional surfaces, and a particle col- 
liding with a wall is analogous to a particle colliding with a particle of infinite mass moving at the 
velocity of the wall. The frictional properties of the walls are modeled using a single parameter, 
the coefficient of tangential restitution for particle-wall collisions (J3 W ), which is defined as the 



3 



fraction of relative tangential momentum transmitted from a particle to the wall during a particle- 
wall collision. The wall-roughness is controlled by choosing a specific value of (3 W : while (3 W = 1 
corresponds to a fully smooth wall, (3 W — corresponds to a fully rough wall for which the dissi- 
pation of energy at walls is maximum and there is no relative tangential slip between the particle 
surface and the wall upon a wall-particle collision. 

Apart from the wall-roughness parameter (3 W , the granular Poiseuille flow is governed by three 
control parameters: the average volume fraction (z/), the coefficient of normal restitution (e„) and 
the channel width (W/d). It may be noted that the gravitational acceleration (g) does not appear 
explicitly as a control parameter since it is used as a reference scale for time (y/W/g), velocity 
(y/Wg) and other mean fields. In the present simulation, we have fixed N = 900 and W/d = 31 
and varied the channel length (L/W) to change the average volume fraction, 

v = vriV/4 {L/W) (W/d) 2 . 

(The robustness of reported results was checked by increasing the number of particles by four- 
fold.) The system is initially allowed to attain a statistically steady state for which the stream- wise 
velocity (JJ X ), volume fraction (v) and granular temperature (T), remain invariant in time, but have 
spatial variations along the wall-normal direction (y). All the statistics presented in this paper are 
computed 'bin-wise' by dividing the channel into 18— bins along the wall normal direction, and 
collecting the data in each bin over about 300000 collisions per particle after reaching the steady 
state. An increase in the number of bins to 31 did not alter the results; a few bins are indicated by 
arrows in Fig. 1(b), with bin = 1, 18 being located adjacent to the walls and bin = 9, 10 at the 
center of the channel. It is to be noted that u x = c x — U x (t) and Uy — Cy tire the particle fluctuating 
velocity in x— and y— directions, respectively, over the instantaneous mean velocity; here U x (t) 
and U x =< U x (t) > are the 'instantaneous' mean and 'time- averaged' mean streamwise velocity, 
respectively. 

Figure HJa) shows the probability distribution function of the fluctuating streamwise velocity 
(u x ) for dilute-to-dense flows (0.015 < v < 0.56) in different bins. The wall-roughness has 
been set to (3 W = 0.9 for smooth walls, and the restitution coefficient to e n = 0.99 for quasi- 
elastic particle collisions. (The distribution function of the fluctuating transverse velocity, u y , 
looks similar, and hence not shown.) The local (bin-wise) distribution functions on only one- 
side of the channel-centerline are presented as the distributions on the other side is the same; 
however, in some cases, the distributions on both sides are presented when they differ. Note that 
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FIG. 1 : (Color online) (a) Distribution function of u x for a range of volume fractions in the quasi-elastic 
limit (e n = 0.99) for smooth ((3 W = 0.9) walls in different bins. The dashed curve indicates a Gaussian. 
Left inset shows the variation of Knudsen number, Kn, with volume fraction, (b) The mean velocity (U x ), 
granular temperature (T) and volume fraction {v) profiles across the width of the channel for e n = 0.99 
and j3 w = 0.9 at different volume fractions. The arrows near the left-ordinate indicate the locations of some 
bins. 

the horizontal axis in the velocity distribution plots is scaled by a t = ■sj (it?), where the index 
i denotes the coordinate direction, and the vertical axis has been scaled such that P(0) = 1. It 
is remarkable that the velocity distribution function in all bins remains a Gaussian for a wide 
range of densities (p < 0.6). This is a surprising result, especially in the dilute limit, since the 
Knudsen number [see left inset in Fig. [TJ a)], which is the ratio between the mean free path and 
the channel width, Kn = lf/W, increases with decreasing v and becomes of 0(1) in the dilute 
limit, signalling the onset of rarefied flow. Even in this rarefied regime, the velocity distribution 
function remains a Gaussian in granular Poiseuille flow with smooth walls. From the profiles of 
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temperature (T), mean velocity (U x ) and volume fraction (u) in Fig.[Qb), we observe that the mean 
field quantities develop considerable gradients along y— direction with increasing density (and this 
is more pronounced for U x , see lower inset in Fig. lb), however, they are almost constant over the 
width of a bin. The mean-field gradients do not seem to play any role in determining the form of 
'local' velocity distribution functions as long as the walls are smooth and the particle collisions 
are quasielastic. 
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FIG. 2: (Color online) (a) Effect of restitution coefficient, e n , on P(u x ) at v = 0.56 for smooth walls 
(J3 W = 0.9). The velocity distribution in bin = 1, upper inset, remains a Gaussian, the distribution in 
bin = 5, lower inset, develops asymmetric tails and the distribution in bin = 9 makes a transition from 
Gaussian to a power-law tails with increasing dissipation. Left Inset: Variation of power-law exponent, a, 
with e n . (b) Same as in panel (a) but for P(u y ) in log-log scale to discern the power-law behviour. 



For a dense flow, v = 0.56, with smooth walls (3 W = 0.9, Fig.[2ta) shows the effect of collisional 
dissipation, e n , on the form of P{u x ) distribution in different bins. The velocity distribution in 
bin = 1 (adjacent to wall), the upper inset in Fig. (2a), remains close to a Gaussian, irrespective of 
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the value of e n . The P{u x ) distribution in bin = 5 (between the wall and the channel-centerline), 
the lower inset in Fig. |2ja), becomes asymmetric with increasing dissipation, with an enhanced 
probability of negative velocities (negative skewness). At the center of the channel (bin = 9), 
the high velocity tails of the distribution function undergo a transition from a Gaussian to a non- 
Gaussian (with over-populated tails) with increasing collisional dissipation, as seen in the main 
panel of Fig. 2(a). The analogue of Fig. 2(a) for the fluctuating transverse velocity u y is shown 
in Fig. Hb). The transition of P(u y ) in the channel-centerline is very similar to that seen in 
P(u x ). The distribution in bin = 1 remains Gaussian and is not shown here, instead we include 
the distribution in bin = 14. Note that P(u y ) distribution in bin = 5 and bin = 14 are mirror 
images since these bins are symmetrically located about the channel-centerline. The appearance 
of asymmetric distribution functions in two shear-layers (with increasing dissipation-level) could 
be a consequence of density- waves Ill2ll in narrow shear-layers. (Within the plug-region, however, 
the local distribution functions are slightly affected by such asymmetries.) This issue is relegated 
to a future study. 

For parameter values as in Fig. 2, the profiles of the mean field quantities are displayed in 
Fig. I^a) which shows the emergence of a 'plug' around the channel-centerline (with negligible 
gradients in U x , T and v) with decreasing e n , and two 'shear-layers' adjacent to two-walls with 
steep gradients in U x , T and v. To pinpoint the role of e n on distribution functions, first we study its 
effect on the pair correlation function and the spatial velocity correlations. The velocity correlation 
function is defined as follows: 

(ui(x)uj(x + 5xj) 
ij ~ MOH(O)} ' 
where the indices i,j denote the coordinate directions. In the quasielastic limit, the pair correla- 
tion function, the lower inset in Fig.[3tb), shows a liquid-like structure in all bins and the velocity 
correlation is close to zero (not shown). The signatures of plug-formation with increasing dis- 
sipation show up in the pair correlation function (the main panel of Fig. |3h) which indicates a 
transition from a liquid to a crystal-like structure in bin = 9 at e n = 0.8. With increasing den- 
sity correlation in bin = 9, the velocity correlation C xx also becomes strong as shown in the 
upper inset of Fig. Ob). (It is interesting to note that the velocity correlation is negative beyond 
a certain correlation length, x/d ~ 10, which is an indicator of circulatory-type motion [13] for 
the fluctuating velocity field.) At e n = 0.8, the pair correlation function outside the plug region 
(bin = 1,5) shows a gaseous structure and the C xx correlation is weak/absent (see upper inset in 
Fig. [2b). Clearly, the enhanced density and velocity correlations around the channel-centerline are 
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responsible for the emergence of non-Gaussian tails with increasing dissipation (Fig. 2). 
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FIG. 3: (Color online) (a) Variations of U x , T and v along y, with parameter values as in Fig. 2. (b) 
Main Panel: Pair correlation function, g{x), in different bins at v = 0.56 for e n = 0.8 with smooth walls 
(f3 w = 0.9). Upper Inset: Streamwise velocity correlation function (C xx ) with parameter values as in main 
panel. Lower Inset: g(x) in different bins for e n = 0.99 which shows a liquid-like structure. 



With reference to dense flows in Fig. 2, the tails of P(ui) in the plug-region can be fitted via a 
power-law of the form P(ui) ~ u^ a \ with a single exponent, a x ~ 7 and a y ~ 5.5, for a range 
of restitution coefficients, see the left insets in Figs. 2(a) and 2(b). The near-constancy of a, for 
e n < 0.85 is due to the fact that the density within the plug saturates to the close-packing limit 
(u c ~ 0.9) for e n < 0.85 and consequently the other hydrodynamic fields also remain invariant 
there with a further decrease in e n . This weak-variation of on on e n is also implicated in its 
variation with average density. A similar power-law behavior (with a ~ 2.9 — 7.4) has recently 
been observed in experiments of gravity-driven channel flow [9]; however, it is difficult to make a 
direct comparison since the experimental geometry is different (with a sieve at the bottom) and the 
flow corresponds to the dense 'quasi-static' regime. This issue can be resolved by probing very 
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FIG.4: (Color online) Effect of mean volume fraction on P(u x ) for a rough- wall (J3 W = 0.1) with e n = 0.99 
in bin = 1. Upper Inset: bin = 9. Lower Inset: Effect of j3 w on P(u x ) at v = 0.015 in bin = 1. 

dense systems which is beyond the scope of the present paper. 

The effect of wall-roughness on P(u x ) is shown in Fig. HI for a rough wall (3 W = 0.1 with 
quasielastic collisions e n = 0.99. The P(u x ) distribution near the wall (bin = 1, the main panel) 
develops a single-peak asymmetric structure at all densities, with its peak being positioned at 
some negative velocity. Near the centerline (bin = 9, upper inset), however, P(u x ) remains 
asymmetric only at low densities and becomes a Gaussian at larger densities. The corresponding 
P(u y ) distribution (not shown) remains a Gaussian at all densities in the quasielastic limit. For 
dense flows with rough walls, both P{u x ) and P(u y ) develop power-law tails with decreasing e n 
around the channel-centerline as in Fig. 2 for smooth walls; the wall-roughness did not influence 
the associated power-law exponents. 

Focusing on the dilute regime (y = 0.015), the effect of (3 W on P(u x ) in bin = 1 is shown 
in the lower inset of Figs. SI It is clear that the asymmetry in P(u x ) diminishes with increasing 
wall smoothness ((3 W ) and becomes a Gaussian when the walls are smooth (f3 w = 0.9). This 
wall-roughness-induced asymmetry in P(u x ) is also reflected in the probability distribution of the 
instantaneous particle streamwise velocity c x (not shown), with a single peak in the low velocity 
region for bin = 1- also, P(c x ) for v = 0.015 approaches a Gaussian with increased wall- 
smoothness. Since the particles loose significant amount of tangential velocity while colliding 
with a rough wall in comparison with their collisions with a smooth wall, a peak near the low 
velocity region is expected for rough walls. The greater the loss of tangential velocity at walls, 
the more is the deviation from a Gaussian, and the related asymmetry in Fig. [4] is, therefore, tied 
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to wall-roughness. On the whole, in dilute flows the effect of wall-roughness is felt on the local 
distribution functions throughout the channel, whereas for dense flows only the near-wall region 
is affected. 

In conclusion, the local velocity distribution functions in a granular Poiseuille flow (GPF) 
with smooth walls remains Gaussian for a wide range of densities for nearly elastic collisions 
(e n — > 1) which, in turn, suggests that the GPF (with smooth walls) could serve as a prototype 
'non-equilibrium steady state' to derive constitutive models starting from the Boltzmann-Engkog 
equation. For dense flows, enhanced density correlations and the related spatial velocity correla- 
tions are responsible for the emergence of power-law tails with increasing collisional dissipation 
around the channel-centerline (irrespective of wall roughness) where the flow undergoes a transi- 
tion from a liquid-like to a crystal-like ('plug') structure in the same limit. For a rough wall, the 
near-wall distribution functions are significantly different from those in the bulk at all densities 
which calls for a relook at the derivation of granular boundary conditions [3D. 
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